Create Map

با استفاده از داده های زلزله ها در ایران و جهان به سوالات زیر پاسخ دهید.

---

کارهای اولیه

library(stringr)
library(readr)
library(dplyr)
library(tidyr)
library(highcharter)
library(ggplot2)
library(ggthemes)
library(ggmap)
library(gganimate)
library(plotly)
library(countrycode)
library(knitr)
theme_set(theme_minimal())

setwd('/Users/aminrakhsha/Documents/University/Term4 - 97spring/Data Analysis/HW/HW11')
sequake = read_rds("data/historical_web_data_26112015.rds")
disaster = read_delim("data/disaster.txt", "\t", escape_double = FALSE, trim_ws = TRUE)
iequake = read_rds("data/iran_earthquake.rds")
earthquakes = read_csv('data/worldwide.csv') %>% filter(type == 'earthquake')

۱. با استفاده از داده های historical_web_data_26112015.rds و استفاده از نمودار پراکنش سه بعدی بسته plotly نمودار طول، عرض و عمق زلزله ها را رسم نمایید. علاوه بر آن بزرگی هر نقطه را برابر بزرگی زمین لرزه قرار دهید.

plot_ly(sequake, x = ~Longitude, y = ~Latitude, z = ~(-Depth), size = ~(10 * Magnitude)) %>%
  add_markers() %>%
  layout(scene = list(xaxis = list(title = 'Longitude'),
                      yaxis = list(title = 'Latitude'),
                      zaxis = list(title = 'Depth')))

---

۲. پویانمایی سونامی های تاریخی را بر حسب شدت بر روی نقشه زمین رسم نمایید.(از داده زلزله های بزرگ استفاده نمایید.)

از پکیج gganimate استفاده میکنیم. نتیجه را در زیر میبینید:

tsunamies = disaster %>% 
  filter(FLAG_TSUNAMI == "Tsu")
# bbox <- c(left = -170, bottom = -60, right = 170, top = 80)
# worldMap = get_stamenmap(bbox, zoom = 3, maptype="toner-lite")
# saveRDS(worldMap, 'worldMap_toner_lite.rds')
# worldMap = readRDS('worldMap_toner_lite.rds')
# p = ggmap(worldMap, extent = "device") +
#   geom_point(data = tsunamies, aes(x = LONGITUDE, y = LATITUDE, frame = YEAR, color = 'red'))
# p = gganimate(p, "q2_1.gif")

worldMap = map_data('world')
p = ggplot() +
  geom_polygon(
    data = worldMap,
    aes(x = long, y = lat, group = group),
    fill = "gray70",
    color = 'gray80'
  ) +
  coord_fixed() +
  geom_point(
    data = tsunamies,
    aes(x = LONGITUDE, y = LATITUDE, frame = YEAR, color = EQ_PRIMARY)
  ) +
  scale_color_gradient2(low = '#59f442', mid = '#f4f142', high = '#f22626', midpoint = 5.5) +
  theme(axis.line=element_blank(),
        axis.text.x=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks=element_blank(),
        axis.title.x=element_blank(),
        axis.title.y=element_blank()) +
  guides(color=guide_legend(title="Magnitude"))
# p
# p = gganimate(p, "q2_2.gif", ani.height = 350, ani.width = 800)

---

۳. نمودار چگالی دو بعدی زلزله های تاریخی ایران را رسم کنید.( از داده iran_earthquake.rds و لایه stat_density_2d استفاده نمایید).

# iranMap = get_map("Iran", zoom = 5)
# saveRDS(iranMap, 'iranMap.rds')
iranMap = readRDS('iranMap.rds')
ggmap(iranMap) +
  stat_density_2d(data = iequake, aes(x = Long, y = Lat, fill = ..level..), geom = "polygon") +
  guides(fill = F)

---

۴. احتمال اینکه در ایران در پنج سال آینده زلزله به بزرگی هفت ریشتر رخ دهد را محاسبه کنید. (از احتمال شرطی استفاده کنید.)

ابتدا با اجتماع گرفتن تمام داده‌ها لیست تا حد ممکن کاملی از زلزله‌های ایران به دست می‌آوریم. تنها زلزله‌های ۱۰۰ سال اخیر که داده‌ی کاملی از آن در دسترس داریم را در نظر می‌گیریم. با استفاده از این لیست توزیع فاصله‌ی میان زلزله‌های بزرگ متوالی به دست می‌آوریم. سپس از روی این توزیع احتمال این که فاصله‌ی زلزله‌ی بعدی با آخرین زلزله کمتر از ۶ سال باشد به شرط این که حداقل یک سال باشد (آنچه تا کنون رخ داده است) را محاسبه می‌کنیم.

iranPowerfullEQs = disaster %>% 
  filter(EQ_PRIMARY >= 7) %>% 
  filter(COUNTRY == 'IRAN') %>% 
  select(year = YEAR, month = MONTH)

iranPowerfullEQs = iequake %>% 
  filter(Mag >= 7) %>% 
  mutate(year = OriginTime %>% as.POSIXlt() %>% .$year + 1900,
         month = OriginTime %>% as.POSIXlt() %>% .$mon + 1) %>% 
  select(year, month) %>% 
  union(iranPowerfullEQs)

iranPowerfullEQs = earthquakes %>% 
  filter(mag >= 7, str_detect(tolower(place), "iran")) %>% 
  mutate(year =  as.POSIXlt(time)$year + 1900,
         month = as.POSIXlt(time)$mon + 1) %>%
  select(year, month) %>% 
  union(iranPowerfullEQs)

iranPowerfullEQs = iranPowerfullEQs %>% filter(year > 1900) %>% arrange(year, month)
diffs = iranPowerfullEQs$year * 12 + iranPowerfullEQs$month
diffs = (diffs - lag(diffs))[-1]

prob = sum(diffs <= 72 & diffs >= 12) / sum(diffs >= 12)
prob

## [1] 0.5714286

---

۵. بر اساس داده های زلزله های بزرگ ابتدا تعداد و متوسط کشته زلزله ها را بر حسب کشور استخراج نمایید. سپس نمودار گرمایی تعداد کشته ها را بر روی کره زمین رسم نمایید.(مانند مثال زیر!)

برای بهتر مشخص شدن تمام کشور‌ها، لگاریتم میانگین کشته‌ها را هم در یک نمودار دیگر رسم می‌کنیم:

disastersByCountry = 
  disaster %>% 
  group_by(COUNTRY) %>% 
  summarise(count = n(), 
            avg.deaths = round(mean(TOTAL_DEATHS, na.rm = T)),
            log.avg.deaths = log10(round(mean(TOTAL_DEATHS, na.rm = T))))
disastersByCountry[tolower(disastersByCountry$COUNTRY) == 'czech republic', 'COUNTRY'] <- 'Czechia'
disastersByCountry[tolower(disastersByCountry$COUNTRY) %in% c('usa', 'usa territory'), 'COUNTRY'] = 'United States'
disastersByCountry[tolower(disastersByCountry$COUNTRY) %in% c('uk', 'uk territory'), 'COUNTRY'] = 'United Kingdom'

codes = codelist %>% 
  mutate(COUNTRY = tolower(country.name.en)) %>% 
  select(code = iso2c, COUNTRY)

disastersByCountry = disastersByCountry %>% 
  mutate(COUNTRY = tolower(COUNTRY)) %>% 
  left_join(codes) %>% 
  drop_na()


hcmap(data = disastersByCountry,
      value = 'avg.deaths',
      joinBy = c('hc-a2', 'code'),
      name = "Average Deaths") %>% 
  hc_title(text = "Average Deaths In Disasters")

hcmap(data = disastersByCountry,
      value = 'log.avg.deaths',
      joinBy = c('hc-a2', 'code'),
      name = "Average Deaths") %>% 
  hc_title(text = "Logarithm of Average Deaths  In Disasters")

---

۶. با استفاده از داده لرزه های بزرگ و به وسیله طول، عرض، شدت، عمق مدلی برای پیش بینی تعداد کشته های زلزله بیابید.

یک مدل خطی تعمیم یافته برای پیشبینی کشته‌ها می‌سازیم. توزیع تعداد کشته‌ها را پواسون در نظر می‌گیریم و از این خانواده استفاده می‌کنیم. نتیجه را در زیر می‌بینید:

# library(h2o)
# h2o.init()  

model = h2o.glm(
  x = c('LONGITUDE', 'LATITUDE', 'FOCAL_DEPTH', 'EQ_PRIMARY'),
  y = 'TOTAL_DEATHS',
  family = 'poisson',
  training_frame = as.h2o(disaster),
  nfolds = 5)

## 
  |                                                                       
  |                                                                 |   0%
  |                                                                       
  |=================================================================| 100%
## 
  |                                                                       
  |                                                                 |   0%
  |                                                                       
  |=================================================================| 100%

summary(model)

## Model Details:
## ==============
## 
## H2ORegressionModel: glm
## Model Key:  GLM_model_R_1531990332305_1 
## GLM Model: summary
##    family link                              regularization
## 1 poisson  log Elastic Net (alpha = 0.5, lambda = 1.8179 )
##   number_of_predictors_total number_of_active_predictors
## 1                          4                           4
##   number_of_iterations training_frame
## 1                    4       disaster
## 
## H2ORegressionMetrics: glm
## ** Reported on training data. **
## 
## MSE:  846409368
## RMSE:  29093.12
## MAE:  7052.057
## RMSLE:  4.681913
## Mean Residual Deviance :  20391.21
## R^2 :  0.007105884
## Null Deviance :40803424
## Null D.o.F. :1592
## Residual Deviance :32483192
## Residual D.o.F. :1588
## AIC :32492210
## 
## 
## 
## H2ORegressionMetrics: glm
## ** Reported on cross-validation data. **
## ** 5-fold cross-validation on training data (Metrics computed for combined holdout predictions) **
## 
## MSE:  876661972
## RMSE:  29608.48
## MAE:  7164.056
## RMSLE:  4.673112
## Mean Residual Deviance :  21456.85
## R^2 :  -0.02838242
## Null Deviance :49727660
## Null D.o.F. :1592
## Residual Deviance :34180762
## Residual D.o.F. :1588
## AIC :34189780
## 
## 
## Cross-Validation Metrics Summary: 
##                               mean           sd   cv_1_valid  cv_2_valid
## mae                      7123.5073    551.82666     6416.725   6273.0815
## mean_residual_deviance   21217.396     4544.921    13384.006   17942.674
## mse                    8.6832653E8 5.32631808E8 3.19777408E8 5.0635072E8
## null_deviance            9945532.0    2991722.8    5158563.0   7118160.0
## r2                     -0.18901128   0.25534597   -0.9085955  0.02662413
## residual_deviance        6836152.5    1661844.9    4189193.8   5275146.0
## rmse                     27378.805      7704.79     17882.32   22502.238
## rmsle                     4.676582  0.061194442    4.7611065   4.7956023
##                         cv_3_valid   cv_4_valid   cv_5_valid
## mae                      7974.5913    6825.1426    8127.9956
## mean_residual_deviance   30539.188    17315.457    26905.656
## mse                    2.3580224E9  4.8355536E8  6.7392678E8
## null_deviance          1.6503202E7    7732185.5  1.3215549E7
## r2                      0.02360808 -0.049356423 -0.037336703
## residual_deviance        9833619.0    5385107.0    9497697.0
## rmse                     48559.473    21989.893      25960.1
## rmsle                    4.6340857     4.622165    4.5699506
## 
## Scoring History: 
##             timestamp   duration iterations negative_log_likelihood
## 1 2018-07-19 16:52:43  0.000 sec          0          40803423.58038
## 2 2018-07-19 16:52:43  0.014 sec          1          34643003.55043
## 3 2018-07-19 16:52:43  0.017 sec          2          32577754.53718
## 4 2018-07-19 16:52:43  0.022 sec          3          32483676.20582
## 5 2018-07-19 16:52:43  0.025 sec          4          32483192.44679
##    objective
## 1 6771.22861
## 2 5751.49936
## 3 5408.92917
## 4 5393.40725
## 5 5393.33726
## 
## Variable Importances: (Extract with `h2o.varimp`) 
## =================================================
## 
## Standardized Coefficient Magnitudes: standardized coefficient magnitudes
##         names coefficients sign
## 1  EQ_PRIMARY     1.085990  POS
## 2    LATITUDE     0.761332  POS
## 3 FOCAL_DEPTH     0.306407  NEG
## 4   LONGITUDE     0.033481  POS

---

۷. با استفاده از داده worldwide.csv به چند سوال زیر پاسخ دهید. تحقیق کنید آیا می توان از پیش لرزه، زلزله اصلی را پیش بینی کرد؟

برای تشخیص پیش‌لرزه و پس‌لرزه‌های یک واقعه به این صورت عمل می‌کنیم: ابتدا با رند کردن طول و عرض جغرافیایی موقعیت مکانی‌را به ناحیه‌ها تبدیل می‌کنیم. سپس در هر ناحیه هر ۵ روز را یک بازه در نظر می‌گیریم. برای هر ناحیه و هربازه‌ی زمانی زلزله‌ها را مربوط به یک واقعه در نظر می‌گیریم. حال برای بررسی پیشبینی پذیری زلزله اصلی فاصله‌ی زمانی اولین لرزه و بزرگترین لرزه برای وقایع مختلف محاسبه می‌کنیم. اگر توزیع این فاصله به گونه‌ای باشد که بتوان با دقت خوبی زمان زلزله اصلی را فهمید یعنی به یک معنی می‌توان پیشبینی را انجام داد. نمودار توزیع در زیر کشیده شده است. می‌بینیم که چنین ویژگی‌ای ندارد و بعد از اولین لرزه، زلزله‌ی اصلی می‌توان در زمان‌های مختلفی بیاید و نمی‌توان زمان آن را پیشبینی و آن روز را تعطیل اعلام کرد. پس این کار ممکن نیست.

منابع دیگر اینترنتی هم این موضوع را تایید می‌کنند که پیشبینی زلزله کار دشواری است و به این آسانی‌ها نیست!

eq.aprox = earthquakes %>% 
  mutate(longitude = round(longitude / 1.5) * 1.5,
         latitude = round(latitude / 1.5) * 1.5,
         day = as.POSIXlt(time)$yday + as.POSIXlt(time)$year * 365) %>% 
  mutate(day = round(day / 5)* 5) %>% 
  arrange(time) %>% 
  group_by(longitude, latitude, day) %>% 
  summarise(firstForeshockTime = first(time),
            primaryTime = time[which.max(mag)],
            count = n()) %>% 
  filter(count >= 10) %>% 
  mutate(diff = primaryTime - firstForeshockTime)


ggplot(eq.aprox, aes(x = diff)) + geom_density() +
  labs(x = "Time Between The First and Strongest Wave (s)", y = "Density")

---

۸. گزاره " آیا شدت زلزله به عمق آن بستگی دارد" را تحقیق کنید؟ (طبیعتا از آزمون فرض باید استفاده کنید.)

از آزمون کوریلیشن با روش spearman استفاده می‌کنیم. می‌بینیم که طبق تست رابطه‌ی مستقیمی بین این دو کمیت وجود دارد!

testData = earthquakes %>% 
  select(depth, mag) %>% 
  drop_na()
cor.test(testData$depth, testData$mag, method = 'spearman')

## 
##  Spearman's rank correlation rho
## 
## data:  testData$depth and testData$mag
## S = 2.8005e+13, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.2265115

---

۹. میانگین سالانه زلزله ها را بر حسب کشور به دست آورید. آیا میتوان دلیلی در تایید یا رد تئوری هارپ ارائه کرد.

این میانگین در avg.eq.year محاسبه شده است. برای بررسی تئوری توطئه‌ی هارپ فرض می‌کنیم که با اعمال این پروژه به طور کلی شدت زلزله‌های یک کشور از سالی به بعد یا در یک سال خاص بیشتر از سال‌های دیگر می‌شود. برای بررسی این پدیده با تست kruskal شدت زلزله‌های سال‌های مختلف را برای هر کشور بررسی می‌کنیم تا ببینیم آیا شدت زلزه‌ها به سال مربوط است یا نه. لیست کشورهایی که در این تست پی ولیو کمی آورده‌‌اند را در زیر می‌بینید. می‌بینیم که حتی خود آمریکا هم در این لیست هست و بقیه‌ی کشورها هم دشمنی خاصی با آمریکا ندارند! پس به نظر می‌آید نباید این تئوری را جدی گرفت.

us.states = state.name %>% paste(collapse = "|") %>% tolower()

earthquakes = earthquakes %>% 
  mutate(year = time %>% as.POSIXlt() %>% .$year + 1900) %>% 
  mutate(country = str_extract(place,',.*$') %>% str_sub(2, -1)) %>% 
  mutate(country = ifelse(str_detect(tolower(country), us.states), 'United States', country))

avg.eq.year = earthquakes %>% 
  group_by(country, year) %>% 
  summarise(avg.mag = mean(mag, na.rm = T))

countries = unique(earthquakes$country)
results = data.frame()
for (cnt in countries) {
  countryData = earthquakes %>%
    filter(str_detect(country,cnt))
  if(length(unique(countryData$year)) < 2)
    next
  testResult = kruskal.test(mag ~ year, data = countryData)
  pvalue = testResult$p.value
  results = rbind(results, data.frame(cnt, pvalue))
}

harped = results %>% filter(pvalue < 0.001) %>% rename(country = cnt)

kable(harped)

country	pvalue
United States	0.0000000
U.S. Virgin Islands	0.0000000
Puerto Rico	0.0000000
Chile	0.0000428
Dominican Republic	0.0000000
Indonesia	0.0000000
Peru	0.0001795
Papua New Guinea	0.0000000
Taiwan	0.0000454
British Virgin Islands	0.0000000
Dominica	0.0000000
Guinea	0.0000000

---

۱۰. سه حقیقت جالب در مورد زلزله بیابید.

به طور کلی زلزله‌ها بیشتر از سونامی‌ها قربانی می‌گیرند.

با تست Wilcoxon تعداد کشته‌های زلزله‌ها و سونامی‌های با شدت حداقل ۶ ریشتر را بررسی کردیم. طبق این تست قربانی‌های زلزله‌ها بیشتر از سونامی‌ها هستند.

testData = disaster %>% 
  filter(EQ_PRIMARY > 6) %>% 
  mutate(FLAG_TSUNAMI = ifelse(is.na(FLAG_TSUNAMI), "Regular", "Tsunami")) %>% 
  mutate(class = floor(EQ_PRIMARY))

regular = testData %>% filter(FLAG_TSUNAMI == "Regular", !is.na(TOTAL_DEATHS)) %>% .$TOTAL_DEATHS
tsunami = testData %>% filter(FLAG_TSUNAMI == "Tsunami", !is.na(TOTAL_DEATHS)) %>% .$TOTAL_DEATHS

wilcox.test(regular, tsunami, alternative = "less")

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  regular and tsunami
## W = 105730, p-value = 0.01501
## alternative hypothesis: true location shift is less than 0

علیرغم رشد جمعیت جهان و شهر‌ها تعداد کشته‌های زلزله‌ها در حدود صد سال اخیر کاهش بافته است.

برای بررسی این موضوع از تست کوریلیشن میان سال و تعداد کشته‌های زلزله‌ها استفاده می‌کنیم. می‌بینیم که رابطه‌ی عکسی میان این دو برقرار است. این درحالی است که با رشد جمعیت جهان و شهرها قاعدتا یاد تعداد کشته‌ها یبشتر شود ولی گویا رشد تکنولوژی ساخت و ساز و روش‌ها مقابله با فاجعه باعث شده است که تعداد کشته‌ها کاهش یابد.

testData = disaster %>% 
  filter(YEAR > 1900)

cor.test(testData$YEAR, testData$TOTAL_DEATHS, method = "spearman", alternative = "less")

## 
##  Spearman's rank correlation rho
## 
## data:  testData$YEAR and testData$TOTAL_DEATHS
## S = 323150000, p-value < 2.2e-16
## alternative hypothesis: true rho is less than 0
## sample estimates:
##        rho 
## -0.3225796

با این که زلزله‌ها بیشتر قربانی می‌گیرند، تعدا خانه‌های تخریب شده‌ی سونامی‌ها و زلزله‌ها تفاوتی معنی‌داری ندارد.

مانند گزاره‌ی اول عمل می‌کنیم ولی این بار نتیجه‌ی تست significant نمی‌شود.

testData = disaster %>% 
  filter(EQ_PRIMARY > 6) %>% 
  mutate(FLAG_TSUNAMI = ifelse(is.na(FLAG_TSUNAMI), "Regular", "Tsunami")) %>% 
  mutate(class = floor(EQ_PRIMARY))

regular = testData %>% filter(FLAG_TSUNAMI == "Regular", !is.na(TOTAL_HOUSES_DAMAGED)) %>% .$TOTAL_HOUSES_DAMAGED
tsunami = testData %>% filter(FLAG_TSUNAMI == "Tsunami", !is.na(TOTAL_HOUSES_DAMAGED)) %>% .$TOTAL_HOUSES_DAMAGED

wilcox.test(regular, tsunami, alternative = "less")

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  regular and tsunami
## W = 2998.5, p-value = 0.2263
## alternative hypothesis: true location shift is less than 0